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Using the dynamical cluster approximation, we calculate the correlation functions associated with 
the nearest neighbor bond operator which measure the z component of the spin exchange in the 
two-dimensional Hubbard model with U equal to the bandwidth. We find that in the pseudogap 
region, the local bond susceptibility diverges at T = 0. This shows the existence of degenerate 
bond spin excitation and implies quantum criticality and bond order formation when long range 
correlations are considered. The strong correlation between excitations on parallel neighboring 
bonds suggests bond singlet dimerization. The suppression of divergence for n <~ 0.78 implies 
that tor these model parameters this is quantum critical point which separates the unconventional 
pseudogap region characterized by bond order from a conventional Fermi liquid. 



Introduction. The low doping pseudogap (PG) region 
of the cuprates has remained an issue of great discussion 
and controversy, with experimental data showing anoma- 
lous behavior such as the suppression of spin excitations 
in the susceptibility, a PG in the single-particle spectra, 
and patterns in the STM spectra, among others [1 0. 
Different investigators argued that the PG is related with 
the settlement of order H 0, S, 0, 0] , where the optimal 
doping is in the proximity of the quantum critical point 
(QCP) associated with this order Previously[l] we 
investigated the staggered flux order in the PG region of 
single-band Hubbard model, as proposed by Chakravarty 
et a/Q. Despite the clear evidence of the PG signature 
in both single particle DOS and two particle magnetic 
spectrum, similar to experimental data in cuprates, we 
found no evidence of staggered flux order. 

Here we investigate a different kind of order associated 
with spin bond correlations. Spin bond order states were 
proposed to take place in the PG region 0, 0, EH U|- 
These bond orders require the investigation of four- 
particle susceptibilities, which is presently very diffi- 
cult to calculate with our method. However, while our 
method does not allow an exhaustive investigation of 
bond order states we find compelling evidence that in the 
PG region the bond magnetic degrees of freedom should 
order. 

Investigating the local bond excitation susceptibility 
with the dynamical cluster approximation (DCA) [12l . 
Il3j |. we find evidence of quantum criticality in the 2D 
Hubbard model. We consider the Coulomb interaction 
U to be equal to the bandwidth W = 8t. The DCA is a 
cluster mean-field theory which maps the original lattice 
model onto a periodic cluster of size N c = L\ embed- 
ded in a self-consistent host. Spatial correlations up to 
a range L c are treated explicitly, while those at longer 
length scales are described at the mean-field level. How- 
ever the correlations in time, essential for local critical- 
ity, are treated explicitly for all cluster sizes. We mea- 
sure the fluctuations associated with the nearest neighbor 
bond operator which measure the z-component of spin 
exchange on the bond. We find that there are degenerate 



bond spin excitations in the doping range 0%— « 22% 
corresponding to the PG region, which results in a diver- 
gent local bond susceptibility at T — 0. This divergence 
is caused by ordering in imaginary time rather than the 
more familiar ordering in space, and associated with the 
settlement of long range order at a general phase tran- 
sition. Nevertheless, in the limit N c — ► oo one should 
expect that long range bond correlations will quench the 
entropy and a transition to a state with long range or- 
der will take place 14| , unless a stronger instability such 
as d-wave pairing occurs first. The DCA method, which 
does not allow spatial ordering on distances larger than 
the cluster size, will fail to capture this transition when 
small clusters are considered. However at temperatures 
larger than the ordering temperature the physics would 
be determined predominantly by the local quantum fluc- 
tuations described with DCA. The divergent behavior of 
bond susceptibility is suppressed for doping > 22% im- 
plying that for these model parameters, 22% doping is 
a QCP which separates the unconventional pseudogap 
region characterized by bond order from a conventional 
Fermi liquid. We also find a strong correlation between 
excitations on parallel neighboring bonds, which suggests 
that the pseudogap region is characterized by bond sin- 
glet dimerization. 

Formalism. To solve the cluster problem we use the 
Hirsch-Fye quantum Monte Carlo (QMC) method 15 1 
which is based on a discrete path integral approximation 
with time step At. Hirsch-Hubbard-Stratonovich (HHS) 
fields are introduced to decouple the interaction[16( 

exp(-ArC/n T n i + ArC/(n T +nj.)/2) = i-Tr^e 2 "^""^ 

(1) 

where an Ising HHS decoupling field a = ±1 is intro- 
duced at each spin-time location on the cluster. This 
transforms the problem of interacting electrons to one 
of non-interacting particles coupled to time-dependent 
fields. The fermionic fields are then integrated out, and 
the integrals over the HHS decoupling fields are per- 
formed with a Monte Carlo algorithm. 

All measurable quantities are completely determined 
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To study bond correlations, we define the bond "y" 
operator at time r as 
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FIG. 1: Phase diagram of the Hubbard model for N c =4 and 
U = W = 8t. T c , Tm and T* are the superconducting, anti- 
ferromagnetic and pseudogap temperatures from In the 
marked region, < S < 0.22 doping, which corresponds to 
the PG regime, a divergent local bond susceptibility is found. 



by the HHS fields and the host Green's function. The 
HHS fields contain all information about correlations 
(spin, pair and charge) in space and time. Moreover, 
Hirsch has shown that spin correlations may be dir ectly 
rewritten in terms of the HHS decoupling fields 1(| 17 1 . 
One may interpret the Ising fields as representing the 
fermion spin variables 



(n(i,r) t -n(i,T);)- (l-e^) 1/2 a(i,r) . (2) 

Note that this is an exact mapping, so that all correla- 
tion functions of the HHS fields are, up to a proportional- 
ity constant, equivalent to the corresponding correlation 
functions of S z (i,r) = ^(n(i, t)j — n.(i,r)jj. 

In the DCA the single-particle and two-particle lat- 
tice response functions are calculated with the Dyson 
equation using the irreducible cluster self-energy and ver- 
tices, respectively. Unlike experiments, where in order to 
search for quantum critcality, a magnetic field is applied 
to suppress the superconductivity, in DCA the super- 
conductivity can be suppressed by imposing a normal 
state host. We then look for divergent susceptibilities in 
the normal state indicating phase boundaries. Lattice 
two-particle susceptibilities indicate transitions to anti- 
ferromagnetic (AF) and d-wave superconducting states 
at finite doping according to the phase diagram shown 
in Fig. [T] However ordering associated with more com- 
plex operators, such as valence bond singletsQ, is far 
more difficult to detect with the DCA, since it involves 
complex equations and up to eight-leg irreducible inter- 
action vertices. More feasible calculations involve the 
corresponding cluster susceptibilities, since they can be 
obtained directly in the QMC process. However, these 
cluster susceptibilities are finite size quantities and can 
diverge only at zero temperature (i.e., infinite imaginary 
time when ordering in time occurs). 



B(i,j;T)=<r(i,T)v(j,T)<xS'(i ) T)S*(j,T) 



(3) 



where "i" and "j" label the position in the cluster. For 
simplicity, we also denote with B nn (B nnn ) the bond op- 
erator when "i" and "j" are (next) nearest neighbor sites. 

In the next section we present results for the correla- 
tion functions: 

Xo(T)= J dT(5B{i;i + x,T)6B{i;i + x,0)) (4) 
X±{T)= fdT{SB(i;i + x,T)SB(i;i + y,0)) (5) 
X|| (T) = J dr(6B(i; i + x, r)SB{i + y; i + x + y, 0)} (6) 



where 



SB(i; i + x, t) = B(i\ i + x,r) - (B nn ) 



(7) 



We also measure \s, the susceptibility associated with 
the operator M 

M{t)= ^£ l (B( i ,i + i;r)+.B( J ,z + y;T)) (8) 
Xs (T)= J dr(S(M(r)6M(0)) . (9) 

These correlation functions describe the response of the 
system to an external field which couples with the bond 
operator B nn . The field acts to modify the z-component 
of the nearest neighbor exchange interaction. Depend- 
ing on its sign it decreases or increases the energy of an 
AF bond and has an opposite effect on a FM bond, xa 
describes the local bond response while an d x\\ the 
correlation between nearest neighbor bonds. x.s is a clus- 
ter quantity incorporating spatial correlations within the 
cluster between bonds. 

Results. We first present calculations on a 2 x 2 cluster, 
the smallest cluster capable of reproducing the generic 
features of the cuprate phase diagram. In the doping re- 
gion relevant for high T c cuprates, the Hubbard model 
shows evidence of short range AF correlations. The ex- 
pectation value of (next) nearest neighbor bond operator 
B nn (B nnn ) is negative (positive) and increases with low- 
ering temperature, as one expects for a system with short 
range AF order. The short range AF order is stronger 
at smaller doping. B nn and B nnn versus temperature at 
different fillings are shown for a N c — 4 cluster in Fig. [2] 
-a) and, respectively, -b). 

In the electron density range 1 > n >~ 0.78, the 
temperature dependence of the local bond susceptibility 
shows the existence of degenerate or almost degenerate 
states with different magnitude of their bond value B nn . 
This doping range roughly corresponds to the pseudogap 
region of an N c = 4 cluster, see Fig.Q] As shown in Fig. [2] 
-c) the extrapolation of our data indicate that xo is di- 
verging when T — > 0. Since the lowest temperature we 
can reach is w O.Oli the apparent divergent xo implies, if 
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FIG. 2: (color online) a) ( b)) Nearest (next-nearest) neighbor 
bond expectation value B nn (B nnn ) versus T for different 
fillings n. Short range AF order is present in the system, c) 
Local bond susceptibility xo versus T. xo show a divergent 
behavior when T — * in the pseudogap region indicating 
critical behavior. 




FIG. 4: (color online) Susceptibility x s {T) at n = 0.88 for 
different cluster sizes. The divergence at T — is more 
pronounced for the large N c = 16 cluster. No divergence 
is present for N c = 2. b) Txs(T) when iV c = 16 for two 
different fillings, c) x\\ ~ X± at n = 0.88 for N c = 16. 
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FIG. 3: (color online) a) Txo(T) show a linear behavior and 
extrapolate to a finite value [1q at T = 0. b)/*o versus filling n. 
c) x\\ (T) increases strongly at low T, presumably diverging. 
However the divergence is weaker than ~ i as Tx\\{T) in b) 
shows. XJ-{T) is small and negative at low T. 



not degenerate states, at least bond excitations with an 
energy much smaller than this scale. 

In the pseudogap region (i.e., 1 > n >~ 0.78) xo di- 
verges as ~ A. This can be seen in Fig. [3] -a) where 
we show (black line) T\o versus T at filling n = 0.88. 
Txo displays a linear behavior and, at T = 0, extrapo- 
lates to a finite value, albeit small, /Iq. The ~ A depen- 
dence of susceptibility is consistent with scenarios which 
assume two degenerate configurations, 1 and 2, with dif- 
ferent bond values such that 2^ = B nn (l) — B nn (2). 
It is instructive to draw an analogy with local spin sus- 
ceptibility of a system with independent local moments. 
A free moment is a doubly degenerate problem where 
the spin can be aligned parallel or antiparallel to a par- 



ticular direction. A perturbing magnetic field lifts the 
degeneracy of the two configuration by an amount pro- 
portional to the moment and the magnetic field. Sim- 
ilarly, in our system the perturbing Hamiltonian acting 
on the bond (ij), H ext — hB(i,j), splits the two config- 
urations with an amount proportional to the field h and 
B(i,j)(l) — B{i,j)(2). Thus, the ~ A dependence of the 
local bond susceptibility suggests the existence of two de- 
generate states with different bond magnitude. Of course 
other scenarios compatible with a ~ A like susceptibility 
cannot be excluded. 

The bond correlations are strongest around n w 0.88 
and weak at small and large doping. This can be seen 
both by inspecting Xo(T) in Fig. [2]-c) and by looking at 
the bond moment /Hq versus filling in Fig. [3]-b). At half 
filling, the bond moment extrapolates to zero, since the 
numerical data does not show evidence of divergent sus- 
ceptibilities in the undoped system. At finite doping ^ 
increases with increasing doping displaying a maximum 
at n « 0.88. ^in) decreases with further doping until it 
vanishes at n rj 0.78. 

At low temperatures we find a strong positive correla- 
tion between nearest neighbor parallel bonds and a small 
negative correlation between nearest neighbor perpendic- 
ular bonds. This is shown in Fig. [3]-c. x\\ 1S increasing 
strongly with lowering T, the numerical data indicating 
even a possible divergence when T — > 0, although weaker 
than ~ A characteristic to local bond fluctuations (see 
Fig. [3] -a). The large value of x\\ shows that increasing or 
reducing the antifcrromagnetism on a bond implies a sim- 
ilar effect on the nearest neighbor parallel bond. Whereas 
the correlation between nearest neighbor perpendicular 
bond fluctuations x_l is much smaller and negative. 

Larger cluster calculations are limited to finite tern- 



4 



peratures due to the minus-sign problem present in the 
Hirsch-Fye algorithm. For example, for the values of 
U/W used here, when N c — 16 the minus sign limits 
Hirsch-Fye QMC calculations to 4tf3 < 36 for fillings in 
the pseudogap region. Down to this temperature the lo- 
cal bond correlation xo (T) looks similar to the one calcu- 
lated for N c = 4, but this temperature is too large for a 
reliable extrapolation to T = 0. However, by increasing 
the cluster size and thus incorporating spatial correla- 
tions at larger length scale one expects a decrease of /io 
or even the disappearance of the divergent behavior due 
to the settlement of bond order. This is similar to the 
disappearance of the T — divergence in the local spin 
susceptibility when going from an atom to finite cluster 
due to the settlement of short range AF order. 

In order to investigate critical fluctuations as a func- 
tion of the cluster size we calculate the cluster bond cor- 
relation Xs (T) defined in Eq. [9l While the zero tempera- 
ture divergence of xo (T) is suppressed in larger clusters 
due to correlations between the bonds, Xs{T) would still 
be divergent provided that short ranged order between 
the bonds emerges. The plots in Fig. [4] -a) of Xs(T) for 
N c = 2, N c = 4 and N c = 16 at n = 0.88 indicate 
that unquenched bond fluctuations persists with increas- 
ing cluster size. For N c = 2 (black), Xs(T), which in 
this case coincides with the local bond xo(T), does not 
show divergent behavior at zero temperature. Xs(T) for 
N c = 4 diverges at T = since it contains the divergent 
Xo(T) term. For N c = 16 and at the accessible tem- 
peratures, Xs(T) increases more strongly with decreasing 
temperature than for the N c = 4 case, thus showing an 
even more divergent behavior. Of course, the large value 
of Xs for N c = 16 show the importance of spatial cor- 
relations between bond fluctuations. As for the N c = 4 
cluster, we find that for N c = 16 the correlations between 
nearest neighbor parallel bonds is more important than 
the correlations between nearest neighbor perpendicular 
bonds, as the difference between x\\ an d X-L plotted in 
FigHJ-c shows. 

Notice that even for the N c = 16 cluster the diver- 
gent behavior of Xs(T) ceases when n <« 0.78 as can 
be seen in of Fig. [4] -b, indicating that n w 0.78 is a 
QCP. The smaller (larger) doping region would presum- 
ably correspond to a state with (without) bond order 
when N c — > oo. 

Discussion Without dismissing other possibilities, we 
note that a scenario where the system forms adjacent 
parallel bond singlets fits very well with our results. For 
instance the divergence of local bond susceptibility xo re- 
quires two degenerate states with different B nn . Suppose 
we measure xo on a bond along x direction. A configura- 
tion with adjacent parallel bond singlets along x, such as 
one in Fig. [5] a), has a B nn = — 1, while a configuration 
with adjacent parallel bond singlets along y, such as one 
in Fig. 03b), has a B nn = 0. If these two configurations 
are degenerate they will yield a divergent xo- Moreover 
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FIG. 5: a) Configuration with bond singlets along x direction. 
The marked bond is a superposition of states with AF aligned 
spins, b) Configuration with bond singlets along y direction. 
The marked bond is a superposition of two states with AF 
aligned spins and of two states with FM aligned spins. The 
bond operator, Eq. |3j measured on the bond along x direction 
(marked bond) takes the value B nn — —1 (B nn — 0) for 
configuration -a (-b). If a) and b) are degenerate, the local 
bond susceptibility will diverge oc = when T — > 0. 



the correlation between parallel bond excitations will be 
also divergent and positive, while the correlation between 
nearest neighbor perpendicular bonds will be negative, 
resembling our findings on the 2x2 cluster. 

Conclusions. The behavior of bond susceptibility in 
the PG region of the 2D Hubbard model calculated with 
DCA shows evidence of quantum criticality and implies 
settlement of bond order. Thus, in the the 2x2 clus- 
ter we find divergent local bond susceptibility at T = 0, 
which implies ordering in the imaginary time due to the 
existence of degenerate bond spin excitations. We find a 
strong correlation between excitations on parallel neigh- 
boring bonds, which suggests that the pseudogap region 
is characterized by bond singlet dimerization. We ar- 
gue that the existence of unquenched local zero energy 
fluctuations for small N c implies long range order in the 
limit A^ c — > oo or the intervention of competing phase 
transition. The suppression of divergence for n <« 0.78 
implies that n w 0.78 is a QCP which separates the un- 
conventional pseudogap region characterized by dimers 
from a conventional Fermi liquid. 
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